// @HEADER
//*********************************************************************//
//  SiYuan: A numerical PDE solver                                     //
//  Copyright (2022) YUAN Xi                                           //
//  This Software is released under the BSD 2-Clause license detailed  //
//  in the file "LICENSE" in the top-level SiYuan directory            //
//*********************************************************************//
// @HEADER

#ifndef _ShallowWater_EquationSet_hpp_
#define _ShallowWater_EquationSet_hpp_

#include <string>

#include "Teuchos_RCP.hpp"
#include "Panzer_EquationSet_DefaultImpl.hpp"
#include "Panzer_Traits.hpp"
#include "Phalanx_FieldManager.hpp"


namespace SiYuan {

  template <typename EvalT>
  class EquationSet_ShallowWater : public EquationSet<EvalT>  {

  public:

    EquationSet_ShallowWater(const Teuchos::RCP<Teuchos::ParameterList>& params,
           const int& default_integration_order,
           const panzer::CellData& cell_data,
           const Teuchos::RCP<panzer::GlobalData>& gd,
           const bool build_transient_support);

    void buildAndRegisterEquationSetEvaluators(PHX::FieldManager<panzer::Traits>& fm,
             const panzer::FieldLibrary& field_library,
             const Teuchos::ParameterList& user_data) const;

  private:
      std::string m_dof_names[3];
      double g_;
      int cellDim;
  };

  template<typename EvalT, typename Traits>
  class Residual_ShallowWater : public panzer::EvaluatorWithBaseImpl<Traits>,
    public PHX::EvaluatorDerived<EvalT, Traits>
  {
    public:
      Residual_ShallowWater(
        const panzer::EvaluatorStyle&   evalStyle,
        const panzer::BasisIRLayout&    basis,
        const panzer::IntegrationRule&  ir,
        const double& g );

      void postRegistrationSetup( typename Traits::SetupData d,
        PHX::FieldManager<Traits>& fm);

      void evaluateFields(typename Traits::EvalData d);

    private:
      /**
       *  \brief The scalar type.
       */
      using ScalarT = typename EvalT::ScalarT;
      const panzer::EvaluatorStyle evalStyle_;

    //  int numDims, numBases, numQP;
      int cellDim;
      double g_;

      /**
       *  \brief depth of water bed
       */
      PHX::MDField<const double, panzer::Cell, panzer::BASIS> bed_depth_; 
      PHX::MDField<double, panzer::Cell, panzer::IP, panzer::Dim> bed_grad_;  // grad of bed depth

      PHX::MDField<const ScalarT, panzer::Cell, panzer::IP> h_;
      PHX::MDField<const ScalarT, panzer::Cell, panzer::IP> hu_;
      PHX::MDField<const ScalarT, panzer::Cell, panzer::IP> hv_;

      PHX::MDField<const ScalarT, panzer::Cell, panzer::IP> dhdt_;
      PHX::MDField<const ScalarT, panzer::Cell, panzer::IP> dhudt_;
      PHX::MDField<const ScalarT, panzer::Cell, panzer::IP> dhvdt_;

      /**
       *  \brief A field to which we'll contribute, or in which we'll store,
       *         the result of computing this integral.
       */
      std::vector< PHX::MDField<ScalarT, panzer::Cell, panzer::BASIS> > fields_;


      /**
       *  \brief The name of the basis we're using.
       */
      std::string basisName_;
      std::size_t basisIndex_;

  }; // end of class Residual_ShallowWater
  
}

#include "ShallowWaterEquation_impl.hpp"

#endif
